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A  SYSTEMATIC  APPROACH  TO  HIGHER-ORDER  PARABOLIC 
PROPAGATION  IN  A  WEAKLY  RANGE-DEPENDENT  DUCT 


1.  INTRODUCTION 

This  work  deals  with  a  continuous  wave  (cw;  i.e.,  monochromatic)  wave  field  and  focuses  on  its 
systematic  separation  in  a  ducting  environment  with  weak  range  dependence  into  components  that 
propagate  in  opposite  directions  along  the  duct  via  a  hierarchy  of  uncoupled  parabolic  equations  (PEs). 
Situations  of  this  general  type  occur  for  electromagnetics  in  the  atmosphere  and  acoustics  in  the  ocean, 
when  very  accurate  modeling  is  required  for  fields  that  propagate  over  long  ranges  through  ducts  whose 
range  dependence  is  weak  enough  that  no  energy  is  backscattered. 

To  highlight  the  essentials  with  a  minimum  of  complexity,  the  simplest  nontrivial  version  of  this 
problem  is  addressed:  waves  of  vertical  (z)  displacement  on  a  membrane  stretched  in  the  horizontal  ( x ,  y) 
plane.  Situations  are  considered  in  which  the  membrane's  density  has  (i)  a  sufficiently  strong  maximum 
in  the  transverse  (y)  direction  to  produce  ducting  in  the  longitudinal  (x)  direction,  but  (ii)  only  a  weak 
dependence  on  the  longitudinal  coordinate.  It  is  shown  that  a  series  of  explicitly  energy-conserving 
“pseudounitary”  transformations  can  asymptotically  decouple  the  Helmholtz  equation  order  by  order  into 
PEs  for  left-  and  right-going  ( -x  and  +x )  propagation.  Parallels  are  drawn  both  to  the  one-dimensional 
limit  of  this  problem  (waves  on  an  inhomogeneous  string  [1,  2]),  in  which  this  splitting  technique 
produces  uncoupled  ordinary  differential  equations  (ODEs),  and  also  to  other  more  complex  “virtually” 
one-dimensional  cases  (e.g.,  in  acoustics  and  elasticity  [3,  4]),  in  which  the  decoupled  fields  propagate  via 
a  hierarchy  of  parabolic  partial  differential  equations  (PDEs)  of  similar  form  to  those  obtained  here. 

It  will  be  shown  that  after  the  mtb  transformation,  the  Helmholtz  equation  assumes  the  form 
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the  a'm>  terms  of  which  contain  a  transverse  differential  operator  e  that  is  assumed  to  be  small  in 

a  specified  sense.  When  the  generator  matrix  Gtm)  is  required,  through  m'h  order,  to  be  diagonal 
( bim)  =  c(m)  =0  so  that  £(m)  and  r<m)  are  decoupled)  as  well  as  balanced  ( a{m)  =  0)  and  symmetrized  ( dim) 
has  the  necessary  symmetry  so  that  left-right  inversion  of  the  problem  simply  swaps  £(m)  and  r(m) ),  this 
splits  the  Helmholtz  field  into  left-  and  right-going  components  that  obey  d£{m)/dx  =  d{m)£(m)  and 
dr{m)/dx  =  -d'"nr1"”  ,  respectively.  Owing  to  the  pseudounitarity  of  the  transformations,  the  original  energy 
is  conserved  at  each  order,  m  and  the  two  counterpropagating  components  separately  transport  their 
shares  of  the  total.  The  form  of  dtm>  that  accomplishes  this  splitting  is  shown  to  be  unique.  Through 
order  m  =  2  ,  d(m)  contains  nothing  but  the  m  leading  terms  of  the  formal  Taylor  series  for  — iV  1  +  2e  ,  with 
dm  being  the  “standard”  PE  operator  [5].  Beginning  at  m  =  3  ,  however,  a  phenomenon  originally  noted 
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by  Wurmser  et  al.  [4]  appears:  d(m>  acquires  additional  mlh-ordcr,  non-Taylor-series  terms  involving  x 
derivatives  of  the  refractive  index. 

2.  BACKGROUND 

This  section  presents  background  information  about  the  membrane  concerning  the  dynamics  in 
general,  its  state-space  description,  and  ducted  propagation. 

2.1  Membrane  Dynamics 

At  rest,  the  membrane  occupies  the  r  =  (x,  y)  plane  at  z  =  0  .  In  the  linear-response  regime,  any 
z  displacement  w(r,t)  evolves  according  to  the  2-D  wave  equation3 

pd2lw-r(yd--dT.)w  =  0,  (2) 

with  phase  speed  c  =  sjr/p  ,  in  which  r  and  p  are  the  surface  tension  and  mass  density.  In  our  case,  t 
will  be  a  constant  but  p  may  depend  continuously11  on  F  .  Energy  conservation  is  embodied  in 
df-p+8te  =  0,  in  which  e  =  \p{d,w)~  +\T(dfw)  is  the  energy  density  and  p  =  -z(dtw)(dTw)  is  the 
energy  flux. 

For  the  cw  case,  with  w(r,t)  =  a  (r)e~iat ,  the  complex  amplitude  obeys  the  2-D  Helmholtz  equation 

[d-r-d?+k2)a  =  0,  (3) 

whose  wave  number  k  =  co/c  is,  in  general,  a  function  of  r  .  Energy  conservation  reduces  to  df  ■  p  =  0 ,  in 
which 


P  =  ~ 


ad~a 


(4) 


is  now  the  time-averaged  energy  flux. 

We  use  the  membrane's  minimum  density  p0  =  min (p)  to  define  the  reference  values 
c0  =  yjt/ p0  =  max(c) ,  k0  =  a>/c0  =  min (k) ,  and  A0  =  2 ^/k0  =  max(2)  in  terms  of  which  the  refractive  index  is 
n  =  c0/c  =  A0/A  =k/k0  =  ^p/p0  >  1 .  For  convenience,  we  also  define  e=  -l)  >  0 ,  a  function  of  x  and y 
that  vanishes  identically  if  the  membrane  is  completely  uniform.  Finally,  we  rescale  r  to  k0r  ,  so  that 
distances  are  measured  out  in  reference  wavelengths;  however,  for  simplicity,  we  continue  to  designate 
this  dimensionless  position  vector  “  F  ”.  Thus,  Eqs.  (3)  and  (4)  become 


“  Subscript  notation  is  used  for  partial  derivative  operators;  e.g.,  djQ  =  SQ/d^  . 

b  In  principle,  a  density  p(x,y)  that  is  discontinuous  across  one  or  more  lines  y  =  h(x)  could  be  handled  by  the 
present  approach,  provided  that  h  varied  negligibly  with  x  over  the  space  of  a  wavelength;  however,  this  would 
only  further  complicate  the  presentation.  The  proper  handling  of  situations  with  significant  h  variability  on  the 
wavelength  scale  is  surprisingly  tricky  [3].  For  simplicity,  density  discontinuities  of  all  types  are  explicitly  avoided 
here. 
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(0_  -0.  +l  +  2e)a  =  0  (5) 

and 

p  =  -ia[a*d-a- ad-a*)  ,  (6) 

in  which  a  =  p0c0  (co/2f  is  a  constant.  We  restrict  our  consideration  to  cases  in  which  p(x,y)  has  a 
pronounced  maximum  in  the  vicinity  of  y  =  0  and  tapers  off  to  its  global  minimum  at  large  |>j .  Thus, 
c(x,  y)  has  a  minimum  near  y-  0  ,  which  we  assume  is  strong  enough  to  form  a  significant  duct  along  the 
x  axis. 


2.2  Splitting  the  Helmholtz  Equation 


By  defining  the  transverse  differential  operator 


s  =  e  +  ±d2 


one  can  write  Eq.  (5)  as 


[^+(l  +  2*)]a  =  0. 


(7) 

(8) 


The  operator  1  +  2c  is  simply  n2{x,y)  +  d2  so  that,  despite  its  altered  appearance,  this  is  still  just  the 
Helmholtz  equation  —  an  elliptic  PDE.  It  has  long  been  recognized  [5]  that,  when  the  environment  has 
no  x  dependence  at  all  (i.e.,  when  n  =  n(y)  so  that  the  operators  c2  and  s  commute),  one  can  formally 
factor  Eq.  (8)  into 

{dx-d)(dx+d)a  =  0,  (9) 

left  right 

in  which  [6] 

d  -  -i\J\  +  2e  (10) 

(Tappert’s  Q  operator  [5])  with  some  suitable  definition  for  the  square  root.  The  left  and  right  factors 
commute  and  should  both  produce  solutions.  However,  it  is  difficult  to  know  how  to  extend  this  insight 
to  x-dependent  environments  in  a  way  that  is  both  mathematically  satisfactory  and  physically  appropriate. 
Indeed, 

The  construction,  approximation,  and  understanding  of  the  square-root  Helmholtz  operator 
have  been  major  problems  in  wave  propagation  and  scattering  for,  at  least,  the  past  fifty 
years  in  areas  as  diverse  as  acoustics,  electromagnetics,  optics,  seismics,  and  pure,  applied, 
and  computational  mathematics.  [7,  p.  1637] 

The  difficulty  is  related  to  the  fact  that  d  is  a  pseudodifferential  operator  [8].  Although  considerable 
progress  has  been  made  towards  mathematically  rigorous  results  (Ref.  7  and  citations  therein),  especially 
for  particular  special  cases  [9],  the  general  theory  is  still  under  development. 
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Rather  than  trying  to  contribute  to  those  mathematical  efforts,  the  present  report  attempts  to  shed 
some  light  on  the  matter  through  another  approach,  namely,  adapting  the  technique  of  Foldy-Wouthuysen 
transformations,  which  has  proven  useful  in  a  similar  context  in  the  quantum  physics  of  atomic  and 
subatomic  systems.  In  the  interest  of  clarity,  two  points  should  be  stressed  at  the  outset.  Firstly,  this  is  a 
method  from  mathematical  physics  that  is,  so  to  speak,  more  physics  than  mathematics.  To  the  best  of  the 
author’s  knowledge,  there  has  never  been  a  thoroughgoing  mathematical  critique  of  the  technique.  Its 
standing  among  physicists  derives  from  its  demonstrated  capacity  to  generate  results  in  agreement  with 
experiment  and  its  success  at  predicting  effects  that  were  subsequently  confirmed  in  the  laboratory. 
Secondly,  though  the  technique  arose  in  the  world  of  quantum  physics  and  is  often  described  using  the 
vocabulary  of  that  field,  there  is  actually  nothing  quantum  mechanical  about  it  at  all.  It  is  simply  a 
physically  motivated  technique  for  splitting  an  elliptical  equation  into  a  pair  of  uncoupled  PEs  under 
certain  conditions. 


The  starting  point  is  the  familiar  observation  [5]  that,  for  small  s ,  the  formal  Taylor  series 
- Jl  +  2s  =  1  +  £■- £-2/2  +  •••  yields  physically  reasonable  results.  In  an x- independent  environment,  when  this 
series  is  truncated  at  1st  order,  the  left  and  right  factors  (dx  +  d)a  =  0  reduce  to  (dx  ±U}  +  e  +  \d2y))a  =  0 . 
With  a  =  l  =  e~"L  and  a  =  r  =  e+,xR ,  respectively,  these  become 


Afl  =  +i 

dx 


1  d1  ^ 

6  H - — 

2  dy  j 

1  d2  ^ 

6  H - — 

2  dy 


L 

R, 


(11) 


the  “standard”  PE  formulation  for  propagation  to  the  left  and  right.  Since  L  and  R  are  slowly  varying 
modulations,  this  outcome  is  widely  held  to  be  valid  to  lowest  order  even  for  cases  in  which  e  has  a 
gradual  x  dependence.  It  seems  intuitively  reasonable  that  retaining  more  terms  in  the  Taylor  series  might 
produce  a  similar  but  more  accurate  result.  Indeed,  this  is  what  happens,  up  to  a  point.  However,  it  turns 
out  that  if  energy  is  to  be  scrupulously  conserved,  additional  non-Taylor-series  terms  must  be  included 
beginning  at  the  third  order.  The  Foldy-Wouthuysen  approach  followed  below  is  a  means  of  generating 
these  necessary  departures  from  the  *Jl  +  2s  Taylor  series  for  environments  in  which  the  x  dependence  is 
gradual  in  a  specific  sense. 


2.3  State-Space  Description 

2.3.1  Initial  Representation 

In  terms  of  the  2x1  “state  vector” 

(12) 

composed  of  the  membrane’s  vertical  displacement  a  and  its  longitudinal  slope  a  =  da/dx,  the  2-D 
Helmholtz  equation  (Eq.  (5))  assumes  the  form 

dxi±  =  Gu_,  (13) 


featuring  the  2x2  “generator”  matrix 
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G  = 


0 


1 


-(1  +  0  0 


=  -egx  +  i(l  +  s)a2 


The  second  form  expresses  G  in  terms  of  the  Pauli  matrices  (Ref.  10,  Ch.  XIII,  Sec.  9), 


(14) 


A  o" 

,  (7,  = 

A  1A 

| ,  (7 2  —  j 

"0  -f 

>  ^3  = 

A  0" 

b 

vl  Oy 

A  0  j 

-A 

(15) 


with  the  unit  matrix  a0  included. 


The  dynamics  of  the  membrane  have  been  expressed  in  terms  of  the  2x1  “state  vector”  (Eq.  (12))  and 
the  2x2  “generator”  matrix  (Eq.  (14)).  Beneath  the  superficial  resemblance  to  the  string  problem  [1,2] 
lies  a  crucial  difference:  G  depends  on  the  transverse  coordinate  v  as  well  as  the  longitudinal  coordinate 
x,  and  the  ordinary  function  e  has  been  replaced  by  the  transverse  differential  operator  s  .  Although  the 
quotes  will  be  dropped  hereafter,  this  distinction  should  be  kept  in  mind. 


The  analysis  could  begin  at  this  point;  however,  since  the  eventual  solution  will  emerge  from  a 
perturbation  process  that  relies  on  s  being  “small”  in  some  sense,  it  will  be  simpler  first  to  shift  to  a 
representation  in  which  the  generator  matrix  becomes  diagonal  whenever  e  vanishes.  This  is 
accomplished  through  the  transformation  Uu  =  u,  in  which 


U  = 


J2 


1  1 


(16) 


is  a  unitary  matrix. 

2.3.2  d  ’Alembert  Representation 

In  this  new  “d’Alembert”  representation,  the  state  vector  is 

1 


u  = 


f a  +  icr 

rt 

Ka-ia j 

AA 

and  the  Elelmholtz  equation  retains  the  canonical  form 

dxu  =  Gu 

in  terms  of  the  transformed  generator 


G  = 


-i(\  +  s)  -is 


=  -i(\  +  s)<7.  +  sg2 


is  i{\  +  s)y 

The  longitudinal  and  transverse  components  of  the  energy  flux  p  are 


and 


px  =  aar^u 


p  =  au  '  P  u  , 


(17) 


(18) 


(19) 


(20) 


(21) 
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in  which  the  f  superscript  indicates  a  complex-conjugated  matrix  transpose,  i.e.,  u  =  (f  ,r‘) ,  and  P  and 
P±  are  the  2x2  matrix  operators 


P»  =  -cr, 


(22) 


and 


PL=i{d\-dy)\(cj0+(jl). 


(23) 


The  f  superscript  in  Eq.  (23)  has  a  different  meaning  that  will  be  addressed  later.  In  fact,  dl  is  a 
transverse  derivative  operator  that  acts  to  the  left;  i.e.,  vo\.  =  dv/dy  . 

Superficially,  Eq.  (18)  resembles  a  state-space  formulation  [11];  however,  it  is  a  PDE.  For  a  genuine 
state-space  description,  it  would  have  to  be  replaced  by  an  ODE.  That  is  exactly  what  happens  when  the 
problem  is  reduced  to  the  “string”  limit  [1,  2]:  the  transverse  (y)  coordinate  disappears  from  the  picture, 
and  the  state  vector  evolves  via  an  ODE  that  is  simply  Eq.  (18)  with  the  operator  e  =  e(x,y)  +  ±d2  replaced 
by  the  function  e(x) .  Then,  G  =  -i(\  +  e(x))cr3  +  e(x)cr2 ,  is  an  ordinary  numerical-valued  generator 
matrix,  and  a  well-defined  propagator  governs  the  system’s  evolution0  in  either  direction  along  the  x  axis 
from  any  initial  point  x0 .  For  the  membrane,  however,  such  a  simple  result  is  not  really  achievable  in  any 
exact  sense  because  the  evolution  equation  always  remains  an  elliptic  PDE  in  x  and  y.  A  state-space 
description  for  the  membrane  can  only  apply  in  an  approximate  sense  and  only  in  environments  in  which 
there  is  no  backscatter  (or  at  most  a  kind  of  “virtual”  backscatter  [3]).  In  the  present  work  the  original 
Helmholtz  equation  which,  being  an  elliptic  PDE,  does  not  lend  itself  to  solution  by  numerical  marching 
methods,  is  “split”  into  a  pair  of  parabolic  PDEs  to  which  marching  methods  can  readily  be  applied. 
Naturally,  this  splitting  cannot  be  achieved  exactly.  How  to  do  it  in  a  systematic,  physically  motivated 
way  is  the  crux  of  the  matter. 

2.3.3  Asymptotic  Ducted  Propagation 

For  the  membrane  problem,  we  will  be  especially  concerned  with  the  total  longitudinal  energy  flux, 

Pfx)  =  \^PAxD’)dy  ■  (24) 

In  Dirac’s  bra/ket  notation  (Ref.  1,  Ch.  1,  Sec.  6),  this  can  be  written  as 

Pw  M  =  a  (u(x,  •)  |  f  |  u(x,  •)) ,  (25) 

in  terms  of  the  inner  product 

(u(x,-)\v(x,-))  =  j  td(x,y)v(x,y)dy  ,  (26) 


c  When  e  vanishes  identically  for  the  string  case,  Eq.  (17)  supports  solutions  consisting  of  two  waveforms  (the 
d’Alembert  modes)  that  propagate  without  distortion  in  the  +x  and  —x  senses.  This  makes  it  a  useful 
representation  in  which  to  begin  analyzing  the  membrane  when  e  is  only  weakly  x  dependent. 
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which  involves  both  matrix  multiplication  and  transverse  integration.  The  operator  8fy ,  as  the  notation 
suggests,  is  the  adjoint  of  dy  relative  to  this  inner  product.  The  state  vector  \u(x,y))  satisfies 
dx\u)  =  G\u) .  Its  dual-space  adjoint  (u(x,y)\  satisfies  dx(u\  =  (w|G+  in  terms  of  the  adjoint  operator 
=  +i{\  +  s t)cr3  +  sfa2 ,  where  =e  +  jd2f .  Note  that  d\=-dy  so  that  (d1^  =d2.  Consequently,  s  is 
Hermitian:  s  =  s' .  The  operators  representing  the  energy  flux  components  are  also  Hermitian  (  P  =  l> 
and  P!  =  P  ). 

Integrating  Eq.  (25)  by  parts  yields 

oo  T=co 

dxP\\(x)  =a^u\x,y){(a0  +  cr1)u(x,y)dy  =  -aut(x,y)P1u(x,y)\  =  -Py(x,y)\\=^  ,  (27) 

— oo  y=— 00 

which  is  a  statement  about  energy  conservation:  any  energy  not  flowing  longitudinally  along  the  duct 
must  be  leaking  away  in  the  transverse  direction.  We  will  restrict  consideration  to  ducted  fields  —  those 
that  propagate  paraxially  without  any  leakage. 


3.  FURTHER  TRANSFORMATIONS 

The  essence  of  the  present  method  is  to  apply  a  series  of  transformations  that  progressively 
diagonalizes  the  generator  while  strictly  conserving  the  longitudinal  energy  flux.  In  Section  3.1,  such 
transformations  are  considered  for  a  general  x-dependent  environment.  In  Section  3.2,  the  analytic  form 
is  given  for  an  x-independent  environment. 

3.1  General  Environment 

The  state-space  description  is  analogous11  to  the  evolution  of  a  two-state  quantum  system  via  the 
Schrodinger  equation,  idt  |w)  =  H\u)  (in  units  for  which  h  =  1 ).  The  analog  Hamiltonian  is  H  =  iG  and  x 

plays  the  role  of  time.  One  may  define  the  “matrix  element”  <p{x )  =  (k(jc,-)  I  <?(*,-)  I  “(*,'))  for  an  arbitrary 
operator  c P(x,y) .  Then  cp  obeys 


dx(p  =  (u\G*0  +  0G  +  d>\u),  (28) 

in  which,  as  always,  the  dot  denotes  a  longitudinal  derivative,  |  =  oc /ox  .  In  general,  in  order  for  q>  to  be 
independent  of  x  we  would  need  to  have  G  d>  +  c PG  +  0  =  0.  When  the  operator  0  lacks  any  intrinsic  x 
dependence,  i.e.,  whencP=0,  that  condition  reduces  to  Gf0  +  0G  =  0  ;  i.e.,  to Hf0  =  0H  .  Any  such 
invariant  is  a  constant  of  the  motion.  (Per  quantum-mechanical  usage,  this  term  is  applied  to  both 
(p  and  0  .)  In  quantum  mechanical  problems,  this  always  applies  to  the  unit  operator  0  =  an  so  that  the 

Hamiltonian  is  Hermitian  (Hf  =  H)  and  the  normalization  of  the  wave  function  is  thus  preserved 
( q>  =  ||w||2  =  constant )  during  the  system’s  evolution  in  time.  Any  transformation  M\ii)  =  \u)  to  another 
representation  is  allowed,  provided  it  is  unitary,  i.e.,  M  '  =  M  so  that  the  Hamiltonian  remains  Hermitian 
after  the  transformation  H t  =  H  . 


d  The  analogy  is  to  the  wave  motion  itself,  not  to  any  quantum  phenomena. 
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For  the  membrane  dynamics,  the  constant  of  the  motion  is  instead  O  =  Pn  and  the  Flamiltonian  is 
consequently  pseudo-Hermitian  [4]:  //V)  =  P H  .  (Note  that  Eq.  (27)  explicitly  restricts  the  pseudo- 
Flermiticity  of//  to  ducted  fields.)  Any  linear  transformation  M  |  u)  =  |  u)  is  allowed  provided  it  is  pseudo¬ 
unitary,  Mf/jj  =  PM  ,  to  ensure  that  p  retains  its  status  as  constant  of  the  motion  in  the  new 
representation.  The  transformed  equation  and  longitudinal  energy  flux  are 

dx\u)  =  G\u)  (29) 

and 

p,i  =a(u\P]{\u)  ,  (30) 

in  terms  of  the  transformed  operators  and 

h = M'p,M  (3i) 


and 


G= A+R , 


(32) 


in  which 


A  =  M'GM 


(33) 


and 


R  =  -M~'M  .  (34) 

Without  any  loss  of  generality,  we  adopt  the  exponential  form  M  =  e~B  ,  which  is  pseudounitary  whenever 
iB  is  pseudo-Flermitian;  i.e.,  whenever 


BTai  +  <j3B  =  0  . 


(35) 


We  proceed  by  seeking  a  form  for  B  that  will  produce  a  diagonal  G . 

3.2  Longitudinally  Invariant  Environment 

It  is  well  known  that  when  the  environment  is  independent  of  x,  this  diagonalization  can  be  done 
exactly  [3].  The  correct  general  form  (produced  by  B  =  \y/ax )  is 


M  =  cosh(ir^)cr0  -sinhif^ier, 


^  cosh(|(//)  -sinh(y(//)^ 

-sinh(Y^)  cosh(|^) 


(36) 


in  which  the  parameter  y/  does  not  depend  on  x.  For  arbitrary  i//  ,  (a)  M  is  independent  of  x  so  that  R 
vanishes  and  the  transformed  generator  is  simply  G  =  A ,  and  (b)  the  back-transformation  to  the  initial 
representation,  |m)  =  TJM\u) ,  is  given  by 
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UM  = 


0~2 


exp(-^) 

-/exp(^) 


exp(-{^)' 

*exp(l^)y 


(37) 


For  the  particular  value, 


<//  =  arctanh 


V  l  +  £ 


(38) 


Eq.  (36)  produces  the  diagonal  generator, 


G  = 


d  0 
0  -d 


=  d  a 3 , 


(39) 


I -  £ 

in  which  d  =  -;'vl  +  2^  .  Since  arctanh  -  =  | log (l  +  2s) ,  this  also  results  in 


( 1  +  ff 


UM  =  -= 

V2 


(l  +  2s)H  (\  +  2s)H 
l  l_ 

-i(\  +  2sy  i(\  +  2sy 


(40) 


which  leads  directly  to  a  membrane  z-displacement  of  the  form  a  =  -j=-(l  +  2 s)  O'  +  r) .  When  only  one  of 
the  two  modes  is  present  —  either  £  or  r  alone  —  its  amplitude  is  related  to  a  by  the  WKB  factor 
i(l  +  2^P. 

Although  this  diagonalization  can  no  longer  be  done  exactly  when  the  environment  is  x-dependent,  it 
can  still  be  done  approximately  if  that  dependence  is  “weak”  in  the  specific  sense  that  d"sldx"  ~  0(e"+1) . 
During  this  process  it  is  essential  to  remember  that  although  the  operator  e  certainly  commutes  with  any 
expression  of  the  form  f{s),  it  does  not  generally  commute  with  other  operators  —  notably  its  own 
longitudinal  derivative,  s  =  ds/dx .  Some  preliminary  results  will  be  obtained  in  the  following  section. 


4.  MATRIX  REPRESENTATIONS 

In  Section  5,  a  series  of  pseudounitary  transformations  of  the  form  M  =  e~B  will  be  devised  to 
progressively  modify  the  generator  G  =  A  +  R  by  preserving  the  diagonality  of  the  A  term  while 
progressively  reducing  the  size  of  the  nondiagonal  term  R  .  The  effects  of  these  transformations  are 
formulated  in  terms  of  the  commutator  [ B ] ,  which  is  defined  by  [  B  \<P  =  \  B.d>\  =  B0-&B  for  any  operator 
0  .  The  whole  development  is  based  on  the  properties  of  a  particular  matrix  representation  for  [5] .  That 
representation  is  introduced  here  in  Section  4,  with  particular  emphasis  on  two  important  situations. 
Cases  in  which  the  [ B ]  representation  is  block-diagonal  are  discussed  in  Section  4.1;  block  off-diagonal 
cases  are  treated  in  Section  4.2. 

With  M  =  e~B  ,  Eq.  (33)  becomes 


A  =  elB]G, 


(41) 


in  which  the  exponential  corresponds  to  the  usual  series 
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=Z^r=1+z+iz2+- 

„=n  n ! 


(42) 


and  Eq.  (34)  takes  the  form 

R  =  E([B])B,  (43) 

in  which  E(z)  =  (ez  -1  )/z  ;  i.e., 

co  n- 1 

£(z)  =  Z— r  =  1+iz+iz2+"'  •  (44) 

tl  n! 

Although  Eq.  (41)  is  well  known  (Ref.  14,  Ch.  3),  Eq.  (43)  appears  to  be  an  original  result.  (The  series  in 
Eq.  (44)  was  not  obtained  by  analytical  proof.  It  was  postulated  by  induction  from  the  first  few  terms  and 
then  verified  for  the  first  1 00  terms  using  the  Maple™  symbolic  computation  environment  to  handle  the 
noncommutative  algebra.) 

Expressing  A  and  R  in  the  forms  shown  in  Eqs.  (41)  and  (43)  suggests  that  it  could  be  productive  to 
treat  state-space  operators  such  as  G,  B ,  A,  and  R  as  vectors  in  a  larger  “superspace”  and  to  recast 
objects  like  [ B ]  in  the  role  of  operators  on  that  space.  We  will  call  such  entities  “supervectors”  and 
“superoperators,”  respectively.  (This  parallels  the  usage  in  statistical  mechanics,  where  the  terms  are 
applied,  respectively,  to  the  density  operator  and  the  Liouvillian  in  the  Liouville  equation  [13].)  Since 
any  state-space  operator  can  be  represented  uniquely  by  a  linear  combination  of  the  Pauli  matrices,  a 
convenient  matrix  representation  for  its  supervector  manifestation  can  be  had  by  simply  stacking  its  Pauli 
coefficients  in  a  column  matrix.  For  example,  state-space  operators  of  the  forms 

G  =  aa0  +  bal  +  ca2  +  dcr3 
B  =  pcr0  +  q<Jl  +  ra2  +  sa} 

(whose  a, ■■■,({  and  coefficients  may  involve  transverse  differential  operators)  can  be  represented 

as  supervectors  by 


G  = 


a 

P 

d 

s 

,  B  = 

.... 

b 

q 

yCJ 

vrv 

(46) 


The  entries  in  the  columns  are  arranged  so  that  those  associated  with  the  diagonal  Pauli  matrices  (  a0,  <r3 ) 
and  those  associated  with  the  off-diagonal  ones  (cr,,  <r2)  are  grouped  separately.  (They  are  separated  by 
dotted  lines  in  Eqs.  (46)  and  (47).)  Corresponding  to  this  representation  for  the  supervector  B,  the  matrix 
representing  the  superoperator  [B]  is 


\p\ 

[7] 

M 

M 

L p] 

-i{r) 

t{q) 

[7] 

i{r} 

\P\ 

-i{s} 

M 

- i{q } 

'M 

[p] 

(47) 
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in  which  the  {•}  notation  indicates  an  anticommutator.  {W}@  =  §¥,<&)  =  W&  +  <PW  .  (The  notation  is  a 
standard  one  in  quantum  optics  [14].) 

The  supervector  representations  shown  in  Eq.  (46)  for  G  and  B  can  be  written  somewhat  more 
compactly  as 


G  = 


B  = 


/  p\ 

Q. 


(48) 


in  terms  of  separate  “diagonal”  (0,  3)  and  “nondiagonal”  (1,2)  submatrices 


This  can  also  be  done  with  the  superoperator  representation  in  Eq.  (47): 


in  which 


[*]  = 


V  a) 

U  5) 


[p] 

W1 

a  = 

f  [ql 

[s] 

\pV 

y-i{r} 

>{q}y 

'[q] 

i{r}  ' 

5  = 

J  P] 

-i{s}' | 

Jr] 

- i{q}j 

[p]  J 

(49) 


(50) 


(51) 


Since  a,---,S  are  2x2  matrices,  they  could  also  be  written  as  linear  combinations  of  the  four  Pauli 
matrices.  The  “diagonal”  ones  in  particular  are  simply 


8  =  a 0[p}  +  a 2{s) . 


(52) 


Thus  far,  we  have  allowed  B  to  have  a  general  form.  As  noted,  however,  Eq.  (35)  must  be  imposed  as 
a  restriction  to  guarantee  the  pseudounitarity  of  e~B .  In  terms  of  the  /;,•••, s  coefficients  in  Eq.  (45),  this 
pseudounitarity  condition  reduces  to  a  requirement  for  q  and  r  to  be  Elermitian  and  p  and  .v  anti- 
El  ermitian;  i.e., 


(53) 

(54) 

In  the  remainder  of  this  section,  two  special  cases  are  treated.  In  the  first  one,  the  matrix 
representation  for  [5]  is  block-diagonal  (i.e.,  a  and  /?  vanish).  In  the  other,  it  is  block  off-diagonal  (i.e.,  y 
and  8  vanish).  In  both  of  them,  it  will  be  convenient  to  take  p  =  0  (hence,  P  =  .?|-)  where  |-}  =  (0,1)T)  so 
that  we  always  have  y  =  a,  [.v  |  and  8  =  a2  { .v  j  . 


and 


P  = 


-P, 


q  =q, 
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4.1  Block-Diagonal  [ B ] 


In  order  to  produce  a  block-diagonal  representation  for  [5] ,  let  q  =  0  and  r  =  0  ,  so  that  both  a  and  /? 
vanish: 


As  a  result, 


\B\  = 


/([*])  = 


7(r)  o  ' 

1°  /(*)J 


(55) 


for  any  analytic  function f.  From  the  simple  forms  y  =  crj.v  [  and  5  =  cr2{s} ,  then6 , 


ym=cr0[s]m  and  Sm  =  (T„{s}"‘  for  even  m 
ym  =  tr,  [i]"'  and  8m  =  <x2  {s}"  for  odd  m 


(56) 


Thus,  by  separating  the  even  and  odd  powers  of  z,  via  / (z)  =  /cven (z)  +  /odd  (z) ,  one  has 

=  { CT0./eve„  ((«])  +  T./<xld  ([■»])  0 

l  0  ^„/evra(M)  +  ^/odd(M) 

In  particular,  for  the  function  exp(z)  in  Eq.  (41),  we  have  expeven(z)  =  cosh(z)  and  expodd(z)  =  sinh(z) ,  so 
that 


,  s  _  |  <T>  cosh  ([s])  +  cr,  sinh  ([s])  0 

^  '  0  <t0  cosh  ({s})  +  o' j  sinh  ({>r}) 


5 


(58) 


and  for  the  function  E(z)  in  Eq.  (43),  we  have  £cven(z)  =  sinh(z)/z  =  C(z)  and 
Eoii(z)  =  (cosh(z)-l)/z  s  S(z) ,  so  that 


E([B])  = 


V0C(M)  +  CTl5([^])  0 

v  0  a0C({s})  +  a2S({s}) 


(59) 


Further,  in  terms  of 


rn 

ro) 

i+H 

7 

’  1  H 

[.)■ 

II 

'o 

(60) 


the  diagonal  and  off-diagonal  blocks  of  the  B  supervector  in  Eq.  (48)  are,  respectively,  P  =  |-).v  and 
Q  =  1 0} ,  so  that 


6  This  follows  from  the  Pauli  matrices’  idempotent  property  cr)  =  a0 ;  i.e.,  from  the  fact  that  they  are  all  square  roots 
of  the  identity  matrix  cr0 . 
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B  = 


|o> 


Consequently,  the  diagonal  and  off-diagonal  blocks  of  the  transformed  generator 


G  = 


are 


(61) 


D  =  (cr0  cosh([.?])  +  cr,  sinh([.v]))Z)  +  (|-)C([1?])  +  |+)6'([.v]))i 


(62) 


and 


N  =  (c70  cosh({.?})  +  cr,  sinh(  {.?}))  IV  .  (63) 

The  only  surviving  requirement  for  pseudounitarity  is  s  =  -/  . 

4.2  Block  Off-Diagonal  [B\ 

To  produce  a  block  off-diagonal  [5]  representation,  let  s  =  0  instead.  Since  we  already  have  p  =  0  , 
that  means  that  y  and  5  vanish,  leaving 


[B]  = 


' 0 
P 


(64) 


as  well  as  P  =  1 0)  and  Q  =  q  |  +)  +  r  |  -) .  Thus,  the  surviving  pseudounitarity  condition  is  as  shown  in  Eq. 
(54).Though  this  case  is  more  complicated  than  the  block-diagonal  one,  the  block  off-diagonal  nature  of 
[5]  is  enough  to  reduce  the  diagonal  and  nondiagonal  parts  of  G  to 


^  =  1 


(apy 


ti  (2 n)\ 


D  +  - 


a 


2/7  +  1 


A^  +  - 


Q 


2«  +  2 


JJ 


(65) 


and 


(Pap 

=o  (2 n)\ 


N  +  - 


2/7  + 


-(PD+Q) 


(66) 


The  matrices  a,  /? make  these  more  difficult  to  deal  with  than  the  matrices  in  Eqs.  (62)  and  (63);  however, 
Eqs.  (65)  and  (66)  do  prove  useful  elsewhere  in  the  report. 

In  the  following  section,  a  series  of  pseudounitary  transformations  will  be  applied  to  the  state-space 
formulation  of  the  problem  for  a  ducting  environment.  Transformations  with  block  off-diagonal  [5] 
representations  will  be  used  to  diagonalize  the  generator.  Those  with  block-diagonal  [ B ]  representations 
will  be  used,  when  required,  to  balance  it.  Two  representational  forms  have  been  devised  for  the 
generator: 
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G  =  ( aa0  +d( r3 )  +  (bcrl  +  ca2 )  = 


f  D ^ 

2nd  form 


the  first  form  corresponding  to  its  role  as  an  operator  on  state-space,  the  second  to  its  role  as  a  vector  in 
the  superspace.  “Diagonalizing”  G  means  making  b  =  c  =  0  (i.e.,  N  =  1 0) ),  and  “balancing”  it  means 

also  making  a  -  0  (i.e.,  D  x  |-) ).  The  result  of  this  effort  will  be  a  systematic  hierarchy  of  increasingly 
higher-order  PEs. 


5,  TRANSFORMING  THE  GENERATOR 

In  this  section,  pseudounitary  transformations  are  applied  to  reexpress  the  state-space  equation,  Eq. 
(18)  in  a  series  of  physically  equivalent  representations.  These  will  be  indexed  by  a  superscript  m ,  with 
Z(m)  denoting  a  quantity  Z  in  the  mlb  representation.  The  fundamental  assumption  is  that  d"s/dx"  ~  s"+l . 
Under  this  assumption,  the  approximation  to  Z(m)  that  contains  only  contributions  through  777th  order  will 
be  denoted  Zw.  (Hence,  Zfm>  =  Z['"]  +  0(sm+I) .)  The  transformations  are  designed  so  that  in  the  777th 

representation,  the  777th-order  approximate  state-space  generator,  G["'1 ,  ends  up  diagonal  and  balanced. 
This  splits  the  7?7th-order  dynamics  into  a  pair  of  uncoupled  PEs  for  left-  and  right-going  wave  propagation 
along  the  duct.  We  will  follow  this  procedure  through  m  =  4  ,  laying  out  the  results  here  and  consigning 
the  full  details  to  the  appendix. 


5.1  Order  m  -  0 

The  process  begins  in  the  d’Alembert  representation  (Section.  2.3.2),  which  we  denote  by  777  =  0 .  The 
equation  of  motion, 


0r|w,O))  =  G<O)|w(O)),  (67) 

has  a  generator  G(0) ,  whose  supervector  representation  has  diagonal  and  nondiagonal  parts 
Dm  =  -7(1  +  c)  |  -)  and  N'i])  =  s  |  -} .  Through  0th  order,  this  generator  is  already  diagonal  and  balanced, 

W<»  =  |o)+ow 

,£><»= '  ’ 

so  that  the  0th-order  wave  field  consists  of  left-  and  right-going  modes  that  propagate  via  a  pair  of  PEs  that 
are  coupled  only  at  1 st  order, 


d£(0) 

dx 

dr(0> 

dx 


=  w[0¥0) 


\  +  0(£), 


=  -dl0V0) 


(69) 


with 


W[01=l. 


(70) 
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5.2  Order  in  — 1 

Next  we  attempt  to  construct  a  pseudounitary  transformation 

Mm\um)  =  \ui0)),  (71) 

to  the  Bremmer  (m  =  1 )  representation  [2],  in  which  the  equation  of  motion, 

dx\um)  =  Gm\um),  (72) 

has  a  generator  Gm ,  whose  lst-order  form  G[1)  is  diagonal  and  balanced.  The  exponential  form 
Mi0)  =  exp(-f?(0))  is  successful  at  this,  provided  that  the  [f?(0)]  representation  is  block  off-diagonal  with 
Qw  ~  s  .  In  summary, 


yields 


B 

q 


(0)  (0) _ 

-  q  ’(7X 

(°)  _ir 
2  b 


Nm  =|0)  +  O(£-2) 

iDm  =  j\  +  2s\-)  +  0{s2). 


(73) 


(74) 


Thus  the  1  st-order  wave  field  consists  of  left-  and  right-going  modes  that  propagate  via  a  pair  of  PEs  that 
are  coupled  only  at  2nd  order, 


5l(1) 

dx 

drm 

dx 


=  +dm£(l) 


\  +  0{s2). 


=  -dmrm 


(75) 


The  operator  on  the  right-hand  side  is  given  by 

idm  =l  +  £. 


(76) 


5.3  Order  m  —  2 

We  continue  the  process,  trying  to  construct  a  pseudounitary  transformation 

M(1)|m(2))  =  |w(1))  (77) 

to  an  m  =  2  representation,  in  which  the  equation  of  motion, 

5,|«l2))  =  G(2)|«(2)),  (78) 

has  a  generator  G(2)  that  is  diagonal  and  balanced  through  2nd  order.  The  usual  form  M(1)  =exp(-8(1)) 
works  again,  provided  the  [5(1)]  representation  is  block  off-diagonal  and  Qm  =  i\cr1  ~  s2 .  In  that  case, 
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B 

q 

r 


(1)  =  qa)ax  ■ 

(!)  _  _I_2 
~  2  6 

<»  __1  6 
~  4* 


(1) 


o-2 


(79) 


produces 


N(1)  =|0)  +  O(£-3) 

/D(2)  =y/i+2s\-)  +  0(£3), 


(80) 


so  that  the  2nd-order  wave  field  consists  of  left-  and  right-going  modes  that  propagate  via  a  pair  of  PEs 


—  =  +dl2¥2) 
dx 


ar(2) 

dx 


=  -dl2]r<2) 


+  0(^3), 


which  are  coupled  only  at  3rd  order.  The  operator  is 


idm  =  l  +  £--£2. 
2 


(81) 


(82) 


5.4  Order  m  -  3 


We  attempt  to  use  the  same  procedure  to  construct  a  pseudounitary  transformation 

M(2)|m(3)}  =  |k(2>)  (83) 

to  an  m  =  3  representation,  in  which  the  equation  of  motion  is 

dx\u™)  =  G0)\u«))  (84) 

and  the  generator’s  3rd-order  approximation  G[3]  is  both  diagonal  and  balanced.  This  time,  however,  the 
familiar  form  M<2)  =  exp(-i?<2))  with  the  block  off-diagonal  [5(2)]  representation  and  Q'~2)  ~  £3 ,  i.e., 


B 

q 

r 


(2) 


=  q 


m. 


s3  £ 


{£,£} 


(2) 


(85) 


does  not  succeed  completely.  Instead,  it  produces  an  intermediate  3rd-order  representation 
|z7<3)^  =  expC8,2,)|M<2)^  ,  in  which  the  equation  of  motion, 


Sv  |m(3)}  =  G(3)  | z/(3> )  , 


(86) 


has  a  generator  G(3)  that,  through  3rd  order,  is  diagonal  but  not  balanced: 
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N(3)  =  |0>  +  0(£-4) 

iD0)  =  >/l  +  2s  \ -)  +  i±[e,  e]  |  +)  +  0(s *)  .  (87) 

v - ; - v - ' 

imbalance 

Balance  can  be  provided,  however,  by  a  second  exponential  transformation  exp(-fi'2))|//(3)^  =  |z7(3)^ ,  in 
which  the  [7?<2)]  representation  is  block-diagonal  with  P,2)  ~  s2 ,  i.e., 


B  =  s  o 
s(2)  =  i\£ . 


(88) 


The  result  is  the  anticipated  equation  of  motion  with  a  generator  Gt3)  that  has 

N0>  =  1 0)  +  C>(£-4 ) 

/Z)(3)  =(j\  +  2s-\s}\-)  +  0(s'i)  .  ^ 

Since  the  3rd-order  approximation  Gt31  is  both  diagonal  and  balanced,  the  3rd-order  wave  field  consists  of 
left-  and  right-going  modes  that  propagate  via  a  pair  of  PEs  coupled  only  at  4th  order, 


8l(3) 

dx 

0r(3) 

dx 


■+dl3]l0) 

:-d[V3) 


}  +  0(£4). 


(90) 


The  operator  is  given  by 


idl3]  =  1  +  £ 


1  2  1  3 

—e  +—£ 

2  2 


(91) 


Note  that  m  =  3  is  the  first  order  at  which  idlm]  contains  a  term  (boxed)  that  does  not  originate  in  the 
Taylor  series  for  *Jl  +  2s .  The  overall  transformation  from  the  m  =  2  representation  is 
Mi2)  =  exp(-fi(2))exp(-_S(2’) . 


5.5  Order  m  =  4 


We  now  attempt  to  extend  the  process  to  4th  order,  seeking  a  pseudounitary  transformation, 

M‘3,|m(4))  =  |w(3)},  (92) 

to  an  m  =  4  representation,  in  which  the  equation  of  motion  is 

dx  |w<4))  =  G(4)  |«<4))  (93) 

and  Gm  is  diagonal  and  balanced.  Once  again,  the  exponential  form  M(3)  =  cxp(-/i'3’)  with  block  off- 
diagonal  [5(3)]  and  Qi3)  ~  £4 ,  i.e., 
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B0)  =  gl3,cr1  +  rf3)a2 

?(3)  =-£*+  i  {£,£}  +  j£2  (94) 

rm  =-%{£2,£}-^£££  +  Is¬ 
is  not  completely  successful.  Instead,  it  produces  an  intermediate  4th-order  representation 
|w<4))  =  exp(R(3))|zf3)^ ,  in  which  the  equation  of  motion, 

|m(4))  =  G<4)  |m(4))  ,  (95) 

has  a  generator  G(4)  that,  through  4th  order,  is  diagonal  but  not  balanced: 

Nw  =  |O)  +  0fr?5) 

i&4>  =  ( -|(^2))|-)-*i[^]|+)  +  0(s5 ) .  (96) 

v - v - ' 

imbalance 

However,  balance  can  again  be  restored  by  a  second  transformation  |m(4)^  =  expfB(3,)|w(4)^ ,  in  which  the 
fB<3’]  representation  is  block-diagonal  and  Pm  ~  £3 ,  i.e.,  and 


R,3)  =  5(3)<T3 

^(3)  =-i-k{£’£}- 


(97) 


The  result  is  the  anticipated  equation  of  motion  with  a  generator  Gt4)  that  has 

N{4)  =|o)  +  0(£-5) 

iD{4>  ={^1\  +  2s -\s +  \s2  +1|{^,£:}j|-)  +  0(s5) . 


Since  Gl4]  is  diagonal  and  balanced,  the  4th-order  wave  field  consists  of  left-  and  right-going  modes  that 
propagate  via  a  pair  of  PEs  coupled  only  at  5th  order, 


d/(4) 

—  =  w[4¥4) 

dx 


dr(4) 

dx 


-d[4V4) 


+  0(*5). 


(99) 


The  operator  is  given  by 


Aj 

-H  | 
+ 
r'co 
~  1 

1 

Co 

+ 

II 

1  .. 

- £ 

5  4 

- £  + 

1  .2  3  .  ... 

—  £  +  —  {£,£} 

2  2 

8 

8 

4  16 

(100) 


At  this  order,  id[m]  has  two  contributions  that  are  not  associated  with  the  Taylor  series  for  Vl+~ 2s  :  the 
boxed  3rd-order  term  from  idm  and  a  new  4th-order  contribution  (double  boxed).  Equation  (100)  agrees 
with  Eq.  (51)  of  Ref.  4  and  Eq.  (4.12)  of  Ref.  3.  The  overall  transformation  from  the  m  =  3  representation 
is  M(})  =exp(-R(3))exp(-R(3)). 
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6.  SUMMARY  AND  DISCUSSION 

Single-frequency  “drumhead”  vibrations  on  an  infinite  membrane  have  been  considered  for  the 
situation  in  which  the  membrane’s  density  p(x,y)  has  only  a  weak  dependence  on  x  (the  longitudinal 
coordinate)  but  a  strong  enough  transverse  maximum  at  y  «  0  to  form  a  duct  along  the  x  axis.  For  this 

case,  strictly  energy-conserving  transformations  M{m)  have  been  devised  that  convert  the  Helmholtz 
equation  (in  dimensionless  scaled  x,  y  coordinates)  into  the  form 

j-\u{'n))  =  Glm]\u(m))  +  0(£m+>)  for  0  <  w  <  4  ,  (101) 

in  which  the  2x2  generator  matrix  Gl",]  —  a  function  of  the  transverse  differential  operator  e  =  e+\82 , 
with  e(x,y)  =  j(n2(x,y)-l^>0  —  is  diagonal.  To  777th  order,  Eq.  (101)  decouples  into  a  pair  of  separate 

PEs  for  longitudinal  propagation  without  backscatter  in  either  direction  along  the  duct.  The  essential 
assumption  is  the  smallness  condition  dms/dxm  ~  em+l . 

The  m th  representation  is  related  to  the  0th  (d’Alembert)  one  through  W(m)  | u<m)j  =  |h(0)^  ,  in  which  Wtm) 

is  the  cumulative  transformation,  W<m)  =  M(0)Mm  .  The  d’Alembert  representation  is  related  to 

the  original  representation  composed  of  the  membrane’s  displacement  a  and  longitudinal  slope  a  by 
(/  |f/01)  =  \u) ,  in  which  LJ  is  a  constant  unitary  matrix,  as  shown  in  Eq.  (16).  The  state  vectors  in  the 

original,  d’Alembert,  and  777th  representations  are:  |m)  = 


,  W,0,)  = 


( 


v7  j 


and 


Am) 


(  f(")  I 


..On) 

V7  J 


To  obtain  the  ducted  wave  field  propagating  in  either  direction  to  order  m  : 

(1)  Specify  displacement-and-slope  initial  values  «(x0,y))  at  some  longitudinal  starting  point  x0 . 

(2)  Convert  these  to  initial  values  in  the  777th  representation  via 

(UWlmjy'  (x0,y)  | u(x0,y))  =  | u(m\x0,yj) .  (102) 

(3)  Solve  the  7??th-order  equation  of  motion  to  the  left  or  right,  i.e.,  either 

r)Pm)(r  vl 

— - — —  =  +d{m]{x,y)ftm\x,y)  for  x<x0  given  l{m)(x0,y)  (103) 

dx 


or 

r!pm)(x  vt 

- — - —  =  -d[m\x,y)rim\x,y)  for  x>x0  given  Pm)(x0,y)  .  (104) 

ox 

(4)  Return  to  the  original  displacement-and-slope  representation 

|m(A>0)  =  UWl'n\x,y)\i^m\x,yy  .  (105) 

The  correct  form  of  the  operator  id[m]  for  m  <  4  is  obtained  by  truncating 
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id  =  yjl  +  ls 

Taylor 


--S  +  —  S2  +  —  {£,s}  +  0(s5) 
8  4  16 

' - V - ' 

3rd  4th 


(106) 


at  order  m.  Besides  the  terms  of  the  formal  Taylor  expansion  for  the  square  root,  this  contains  the 
additional  3rd-  and  4th-order  terms  that  are  the  principal  result  of  this  work.  The  effects  of  the  x 
dependence  of  p(x,y)  accumulate  during  step  (3),  and  for  accuracy  in  higher-order  propagation  over  long 
ranges,  the  non-Taylor-series  terms  are  necessary. 

On  the  other  hand,  the  effects  of  the  end-point  transformations  (from  the  initial  representation  and 
back  to  it  in  steps  2  and  4)  do  not  accumulate.  With  that  being  the  case,  a  reasonable  approach  for 
computing  UW[m,(x,y)  is  to  assume  local  x-invariance  at  the  end  points  of  the  propagation;  i.e.,  to 
calculate  it  as  if  n  =  n(y) ,  as  was  actually  the  case  in  Section  3.2.  This  simply  means  using  the  WKB  form 

Eq.  (40)  and  obtaining  ( UW[m] )  '  (x0,y)  and  UWlm](x,y) ,  by  truncating  the  Taylor  expansions  of 


^'=T2 


V2 


(  1  o 


(l  +  2ey 

i(  1  +  2s)“4 

(l  +  2^)4  - 

-i(l  +  2fp 

(1  +  2 s)~* 

(1  +  2  £p 

-i(\  +  2ey 

/(l  +  2p 

(107) 


at  the  appropriate  order. 

This  overall  procedure  rests  on  two  basic  ideas: 

(1)  using  the  Foldy-Wouthuysen  method  of  strictly  energy-conserving  transformations  with  the 
assumption  dnxs  ~  s"+]  to  obtain  the  unique  form  of  the  /r/'-ordcr  decoupled  PEs,  and 

(2)  taking  the  environment  to  be  locally  longitudinally  invariant  at  the  beginning  and  end  of  the 
propagation  so  as  to  exploit  the  WKB  form  for  the  endpoint  transformations. 

The  first  was  conceived  by  Wurmser  et  al.  [4]  who  implemented  it  in  an  ocean  acoustics  application 
involving  a  nearly  range-invariant,  ducting  environment.  The  second  seems  to  have  originated  with  the 
late  F.  D.  Tappert  [15,  5]  who  suggested  that  propagation  in  such  situations  be  done  by  first  extracting  a 
WKB  factor  at  the  source,  then  propagating  the  remaining  auxiliary  field  (the  carrier  of  the  down-range 
energy  flux)  to  the  receiver,  and  finally  including  the  WKB  factor  for  the  receiver.  The  results  obtained 
here  agree  with  corresponding  outcomes  in  Refs.  3  and  4,  where  such  a  procedure  is  followed  in  different 
physical  contexts. 

In  general,  a  numerical  algorithm  that  models  the  higher-order  ( m  >  2 )  parabolic  propagation 
equation  under  the  assumption  that  idl",]  is  simply  the  truncated  Taylor  series  for  J\  +  2c  ,  will  need 
reprogramming  in  order  to  incorporate  the  effects  of  the  non-Taylor-series  terms  in  Eq.  (106).  The 
exception  is  order  m  =  3,  where  the  first  and  largest  of  these  terms  appears.  At  this  order,  Eq.  (106) 
yields  id[i]  =  (\  +  s-\s2  +  p3)-}£  .  Flowever,  since  d"s/dx"  is  simply  d"e/dx"  ,  the  final  term  reduces  to 
-je  ,  which  is  not  a  transverse  differential  operator  but  an  ordinary  function  of  x  and  y.  The  3rd-order 
operator  is  thus  idm  =  (l  +  s -\s2  +  js3 ) - {e  .  Since  s  =  e  +  {d2 ,  the  same  result  can  be  produced  from  the 
Taylor-series  form  idl3]  =  l  +  s-\s2  +\s2  by  making  the  replacement  e->e-je  (and  dropping  4th-  and 
higher-order  terms).  This  can  be  implemented  by  simply  reusing  the  Taylor-series  numerical  algorithm 
with  the  actual  density  p(x,y )  replaced  by  p(x,y)-\d2 p(x,y)/dx2  . 
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This  work  has  presented  a  detailed  analysis  of  one-way  ducted  wave  propagation  for  the  specific  case 
of  a  vibrating  membrane.  The  objective  in  writing  it  has  not  been  to  satisfy  any  perceived  groundswell  of 
community  interest  in  vibrating  membranes  (even  the  author  could  not  deceive  himself  that  far)  but  rather 
to  present  a  full  account  of  a  widely  applicable  solution  method  in  an  uncluttered  context,  with  no  more 
complexity  than  absolutely  necessary,  so  that  the  essentials  would  stand  out  vividly.  This  method  is 
applicable  in  atmospheric  electromagnetics,  optics,  ocean  acoustics,  and  other  wave-motion  disciplines,  in 
which  very  accurate  modeling  is  required  for  fields  that  propagate  over  long  ranges  through  ducts  whose 
range  dependence  is  weak  enough  that  no  energy  is  backscattered. 
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Appendix 


DETAILS  OF  THE  TRANSFORMATIONS 


Al.  INTRODUCTION 

In  this  appendix,  details  of  the  derivations  of  the  transformation  operators  Mlm)  cited  in  the  body  of 
the  article  (Eqs.  (73),  (79),  (85),  (88),  (94),  and  (97))  are  presented.  An  expanded  notation  is  used.  For 
any  quantity  Z  that  can  be  expressed  as  Z  =  Z0  +  Z,  +  Z2  +  •  •  • ,  with  Zk  -  sk ,  the  following  type  of  notation 
is  applied  to  various  portions  of  the  series: 


Z  —  ZQ  +  Z^  +  Z  2  +  Z 3  +  • 


The  subscripts  on  the  terms  indicate  the  ranges  of  ^-orders  that  they  contain.  It  should  be  carefully  noted 
that  these  order  subscripts  are  “applied  last.”  For  example,  the  symbol  7'!'"’  represents  ((p1"0}  the 

order-f  term  in  the  series  for  <b'""  —  not  d&ff  /dx  .  Likewise,  Z(2)  means  (z2)  =  Z0Z2  +Zf  +Z2Z0  ~  s 2 

—  not  (z(2) )  =  Z\  ~  e* . 

For  the  mth  step,  the  objective  is  to  find  a  pseudounitary  transformation  MUn)  that  will  diagonalize  and 
balance  the  transformed  generator  Gtm)  =  Dym)  +  N<m)  through  order  sm  .  In  the  above  notation,  this  means 
a  transformation  that  produces  N$.m)  =  1 0)  and  D$lm)  oc  |  -) .  It  will  turn  out  that  not  only  do  such 
transformations  exist  at  each  order,  but  they  are  actually  unique. 

Two  types  of  pseudounitary  transformations  are  applied  to  carry  the  generator  from  one 
representation  to  the  next.  With  a  block  off-diagonal  transformation,  G  arises  from  G  through  Eqs.  (65) 
and  (66),  which  intermix  G’s  diagonal  and  off-diagonal  components.  With  a  block-diagonal 
transformation,  the  connection  is  provided  by  Eqs.  (62)  and  (63),  which  keep  diagonal  and  off-diagonal 
components  separate.  The  Taylor  Series  for  the  four  functions  that  appear  in  Eqs.  (62)  and  (63)  are 

cosh(z)  =  l  +  |z2  +jjz4  +  ••• 
sinh(z)  =  z  +  jjz*  +jjz5  H — 

C(z)  =  l  +  ^z2+^z4+- 
S(z)  =  jrz  +  ^z3  +fz5  4  . 

After  each  transformation,  it  will  be  found  that  the  dynamics  follow  the  canonical  decoupled  form 
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di{x,y) 

dx 

dr(x,y) 

dx 


= -d(x,y)t(x,y) 
=  +d(x,y)r(x,y ) , 


to  the  appropriate  order  of  approximation,  in  which  d(x,y)  is  the  sum  of  two  parts.  The  first  is  simply  a 
polynomial  function  T(s(x,y )) ,  the  truncated  formal  Taylor  series  for  y]l  +  2 s(x,y) .  The  second  part  is 
the  non-Taylor-series  terms.  Their  details  emerge  from  the  following  paragraphs;  however,  there  is  a 
symmetry-related  restriction  on  their  structure  that  is  clear  at  this  point.  A  typical  form  for  such  a  non- 
Taylor-series  term  is 


sJ(x,y) 


rdFS(x  ,y)'M 
v  dxK  j 


£N(x,y), 


in  which  J,K,L,M,N  are  nonnegative  integers.  Thus, 


8i(x,y) 

dx 


=  ~ \  T ( e(x,y ))  +  AeJ (x, y)  [dKsL (x, y)/ dxK  )“  sN (x,y)  l t{x, y) 


dr(x,y) 

dx 


typical 


same 


}  r(x,y). 


where  A  is  some  constant.  For  simplicity,  this  includes  only  a  single  typical  non-Taylor-series  term.  It 
would  be  trivial  to  extend  the  present  argument  to  the  general  case.  If  the  environment  is  reflected  in  the 
longitudinal  direction  through  e(x,y)  h>  e(-x,  y) ,  this  becomes 


di(x,y) 

dx 

dr(x,y) 

dx 


=  -  \  T (s(~x,  y))  +  As’  (- x,y )  {dKsL  (-x,  y)/ dxK sN  (-x,  y) 


typical 


=  + 


same 


And  if  the  direction  of  the  x  axis  is  reversed  (xh  -jc  ),  this  becomes 


}r(x,y) 


dl(~x,y) 

dx 

dr(-x,y ) 
dx 


=  +  \  T (e(x,  y))  +  A(-  1)KM  sJ  (x,  y)  [dKsL  (x,  y)/ dxK  )“  sN  (x,  y) \l(-x,  y) 


typical 


same 


}r(~x,y)  . 


Left-right  symmetry  demands  that  2(-x,y)  =  r{x,y)  and  r{-x,y)  =  i(x,y) .  For  this,  the  s  exponents 
./. L, N  are  immaterial,  but  it  is  essential  that  the  product  KxM  be  an  even  integer.  Thus,  for  example, 
the  non-Taylor-series  part  might  include  such  terms  as  d2sldx 2 ,  ( ds/dx )2 ,  and  sd2e/dx 2 ,  in  which 
KxM  =  2 ;  however,  it  could  not  contain  terms  like  de2/dx  or  ds2/dx  because  they  have  KxM  =  1 . 

With  this  symmetry  in  mind,  we  proceed  to  derive  the  transformation  operators  M'""  for  m  =  0,1, 2,3  . 
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A2.  ORDER  m  =  0 

The  process  begins  in  the  d’Alembert  representation  as  shown  in  Eq.  (19),  where  the  superoperator 
representation  of  the  generator  G<0)  is  composed  of  D<0)  =  -i(\  +  &■)  |  -)  and  Nm  =  f|-) .  Thus, 


and 


in  which 


D<0)  =  £$  +  D<°> 


N(0)  =  N$ , 


Dm  = 

K?  =  --!-) 


(A.1) 


(A.2) 


(A.  3) 


A3.  ORDER  m  =  1 


For  the  first  step,  to  the  Bremmer  representation  [1],  we  look  for  an  operator  Q(0)  ~  s  to  produce  a 
transformation  Q""  :  G(0)  — >  Gm  that  achieves  the  desired  degree  of  diagonalization,  ^..i,=|0).  Since 
Qm  =0(<i°)),  it  follows  that  «(0,  =«((1°)) ,  /?t0)  =  ,  and  Qm  Thus,  Eqs.  (65)  and  (66)  with  G  =  G(0> 

and  G  =  Gm  yield 

o”..,  -  o;.">„+<W'+ieS) 

(A-4) 


and 


V(1> 

iV(0---4) 


(A.5) 


The  subscripts  facilitate  grouping  the  terms  by  order: 


Dm  =  D<0) 

^(O— 1)  ^(O— 1) 


CD  (2< 

(01A 


(i) 


^u)  =i  (q) 

tV(3)  0 


(°)  +  «(0)D(0) 


(2) 


D(1)  =  1« 

(4)  6  “(1) 


(0)  (0) 
(1)  Md) 


'(1)  (1) 
r(0) 


) 


(JVd)  +  4  PwHo) ) 


(A-6) 


and 
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*SH°> 

V(1)  RW)Om  +O(0) 

ly  (2)  /^(l)  (1)  2^(2) 

Ard)  i  B{0)ai0)  i Nw  +  -  Bw D<0) ) 

iV(3)  2^(l)U(l)  gV(l)  T3Al)  (0)  ) 

Nm  =  1  Bwaw  ( Bm D,0)  +  Ow ) 

JV(4)  6A’(l)a(l)  V-^fl)  (1)  ^(2)  )' 

The  requirement  for  diagonalization  through  order  e  is  that  /V,1,1’  =  1 0) ;  i.e., 


'0" 

+ 

*>(0)}  ^ 

f°)-l 

f°) 

Ay 

vk<0’] 

w 

(A-7) 


(A- 8) 


This  is  equivalent  to  {<7<0)}1  =  £(  and  {r(°’}l  =  0.  Since  {^}1  =  2^  for  any  <j> ,  the  outcome  is  simply 


„«>) 


=  0 


(A.9) 


Because  s  =  for  ducted  fields,  the  pseudounitarity  conditions  in  Eq.  (54)  are  satisfied  at  this  stage. 
Thus, 


g(0,= 


2  b 


i<T  Nm 

2U2iV(l)  * 


(A.  10) 


As  a  result, 


g,0,= 


2  b 


a,0,= 


0 

H£} 


Pw  = 


r\[e\ 


0 

~U£}y 


(A.11) 


(This  is  where  s  first  appears.)  As  expected,  g“”  ~  s2 ,  while  a'y\p{"'  ~e\  i.e.,  g""  =  ,  a'°'  =  , 

and  J3(0)  =  p\[ 2).  Equations  (A.  10)  and  (A.ll)  provide  matrix  representations  for  Eqs.  (A.6)  and  (A.7). 
Those  terms  not  already  obtained  are 


Da)  = 
1J(2) 


A(1)  = 
iv(2) 


f  0  1 

i  „2 

/ 

D(1)  = 

’  ^(3) 

V 

V2^  J 

/  1  •  \ 

/ 

2^ 

„2 

\~£  J 

ItM] 

3 


D(1)  = 

^(4) 


AT(  1)  _  I  48 
JV(4) 


_ 1£ 4 

v  66  y 


(A.  12) 


A4.  ORDER  »i  =  2 


For  the  second  step,  we  look  for  an  operator  g(1)  ~  e2  to  produce  a  transformation  g(1)  :  G(1)  ->  G(2) 
that  diagonalizes  the  generator  through  2nd  order,  N{2)..2)  =  |0) .  Since  g‘"  =  g^J ,  it  follows  that  a'"  =  a®  , 
pm  =  ,  and  g(1)  =  g^ .  Thus,  Eqs.  (65)  and  (66),  with  G  =  G(1)  and  G  =  G(2) ,  yield 


£>(2)  =£>(')  +  a(1)  A211’  +  1«t”  B(l)  Dm 

a^(0---4)  -t^(0---4)  ^  tx,(2)JV(2)  T  2  tx(2)/y(2)x^(0) 


at(2)  _  BmD®  +A(1)  +  6 

J  v  (0---4)  R(2)±j(0--2)  t  iv(2'"41  T  >£( 


(!) 

“(3) 


(A.13) 
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Consequently, 


and 


Ao.),,=A(o)..n=Ao.,,) 

^.3)=^-3,  (A'14) 


<.,)=<i,=l°) 

N(2)  =  Nm  +  6mDm  +Om 

iV(3)  iV  (3)  ^  ^(2)^(1)  ^  ^(3) 
Af(2)  _  at(1)  i  /?0) 


(A- 15) 


The  requirement  for  diagonalization  through  order  s 2  is  that  =  1 0) .  In  the  matrix  representation,  this 
means 


[<7U'] 

[rm] 


i{rm} 

-i{qm} 


20 

y-i 


and  the  solution  is 


qm=~\s\ 


rm=-\s. 


(A- 16) 


(A.  17) 


(Since  s  =  e  is  just  an  ordinary  function,  the  pseudounitarity  conditions  are  satisfied.)  Thus, 


Furthermore, 


0l,,  =  -K+  K*  -  = 


2  U2iV(2)  ' 


(A.  18) 


f?(1)  = 


V  ~7*  V 


«(1)  = 


M[^2] 


,  r  = 


-{[A]  |{A} 


(A.  19) 


Equations  (A.18)  and  (A.19)  produce  the  matrix  representatives  for  the  quantities  in  Eqs.  (A.14)  and 
(A.  15).  Those  not  already  obtained  are 


N{2)  = 

JV(3) 

N<2)  = 

JV(4) 


~\[e\e\ 

ii^4  +  ‘-£2 
V  8  *  T  8  b  J 

4J  _1  •• 

V3  C  4  *  J 


±[s,[s,s}\  +  \{s2,s} 


-f*4 


(A.20) 
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A.5  ORDER  m  =  3 


For  the  third  step,  we  look  for  an  operator  Q<2)  ~  F  to  produce  a  transformation  Q,2> :  G<2>  ->  G(3)  that 
diagonalizes  the  generator  through  3rd  order,  =  |o).  Since  Qm  =  Q‘2’ ,  it  follows  that  a(2>  =  a'3' , 

/?,2)  =  /?(‘32) ,  and  Q(2)  =  0$  .  Thus,  Eqs.  (65)  and  (66)  with  G  =  G(2)  and  G  =  G,3)  yield 


Z)<3)  =  D{1) 

■^(0—4)  ^(0—4) 

)&(3)  =  Af(2)  +  /?(2)  D(2)  +  0(2) 

iV(0-4)  iV(3-4)  ^  /^(3)  -^(0-1)  ^  ^(4)  * 


(A.21) 


Gathering  terms  of  the  same  order,  one  finds  that  all  the  “D”  terms  are  the  same  as  for  m  =  2  while  the 
“iV”  terms  are 


4%=<)-2)=|0) 

Ng>  =  N*>+/r<»D<*> 


(A.22) 


The  requirement  for  diagonalization  through  order  A  is  that  =  |0) .  From  the  matrix  representation, 
one  obtains 


c/2)  =  ff3-{A  r,2)  =\{s,s}=\ds2/dx  . 


(This  satisfies  the  pseudounitarity  condition.)  Thus, 


e(2>=(|,3-^)l+)+l{,,^}|-)= 


^2J  _lc\ 
3  b  8  b 


=  i/j  w(2) 

2  U2iV(3)  * 


Consequently, 


0(2,= 


a<2)  = 


/?(2>  = 


de3 /dx  -\d3sl  dx3  ^ 
v  I  d2s2/dx2 
{  [j£3-j£]  \[d£2/dx\ 

-i\{ds2  /dx}  i{f£3-±s} 

i3j{ds2/dx} 


(A.23) 


(A.24) 


(A.25) 


\[ds2/dx]  -i{\s3-\s}j 

These  produce  the  matrix  representations  for  the  terms  of  Eq.  (A.22).  The  only  one  lacking  is 


Nm  = 

JV(4) 


■jg  [ e ,  [£,  £]]  +  }  {£2 ,  £}  +  { {£,  d£2/dx}  +  \d£3  /dx  -  dx3  ^ 

-2  £4  +j{£,£}+jd2£2/dx2 


(A.26) 


In  summary, 


D{3>  = 

^(0—4) 


l[£,£]-\[£2,£] 

i((VlT 2e) 


'  (0  4) 


(A.27) 
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Although  the  generator  is  diagonal  through  3rd  order,  ( N$..3)  =  |  0) ),  this  time  it  is  not  balanced: 


Z)t3) 

^(0-3) 


ItM] 

-i^l+  2e  j 


(0-3)7 


(A.28) 


Balance  is  achieved  with  a  block-diagonal  transformation,  Pl2) :  G(3'  — >  G(3) .  The  only  workable  option  is 
s  ~  s2 .  (It  is  a  straightforward  though  tedious  matter  to  check  that  all  other  options,  s  ~  s"  with  n*  2 
fail.)  With  that,  Eqs.  (62)  and  (63)  yield 


fl(0--4)  —  fl(3  -4)  +  -2) 

^(0-4)  =  ^(0-4)  +  S 


(A.29) 


To  produce  a$.. 3)  =  0  ,  we  need  a$  +[,v|(/I'3’  =  0  ;  i.e.,  j|>,£]  +  [■?, -ie~\  =  0  .  From  the  general  2nd-order  form 
s  =  0s2  +  <ps  ,  we  find  that  cp  =  ij 8 .  Flowever,  0  is  not  determined  by  the  balance  requirement  (because  s 2 
commutes  with  f).  Thus,  s  =  0s2+i\s,  which  yields  a((3)..4)  = -1|[£2,£-]  and 

4%  =-^Vl  +  2fj  +  /(-A  +  Ode2 jdx+i\s  .  Since  the  left-right  symmetry  requirement  excludes  terms 
of  the  form  ds2/dx  ,  we  must  have  0  =  0.  In  summary, 


N<3)  = 

(0--4) 


-kl£f  £  +  1  {^2  }£  +  f  M  Ss2/dx  - 


\ds 2  /  dx  -  \d3£/dx3  ^ 


- 2s 4  + 1  {s}s  +  |<3  V 2  j  dx 2 


(A.30) 


and 


n(3)  = 

-^(O— 4) 


/ 

—  3  1 

16  1 

'sjl  +  2e) 

V 

K  /( 0- 

(A.31) 


Through  3rd  order,  the  transformed  generator  remains  diagonal  ( IV®. 3)  =  1 0) )  and  is  now  also  balanced: 


Ao-3)  - |-t(Vl  +  2£-)(0  ^  (A.32) 

though,  for  the  first  time,  contributions  appear  in  it  that  do  not  belong  to  the  Taylor  series  for  ~J\  +  2s  . 

A6.  ORDER  m  —  4 


The  fourth  step  begins  with  a  block  off-diagonal  transformation,  Q'3' :  G{3)  ->  G(4) ,  in  which  it  is 
assumed  that  0t3)  -  A .  Since  Q<3)  =  g)3] ,  we  have  a(3)  =  a23’ ,  /?(3)  =  /?(t3' ,  and  Q(3)  =  .  With  m  =  3  in 

Eqs.  (65)  and  (66)  (with  G  =  G(3)  and  G  =  G(4) ),  the  n>  0  terms  and  those  containing  Q(3)  contribute 
nothing  through  4th  order,  so  the  result  is  simply 

Z)(4)  =  p>(3) 

-^(0—4)  -Ly(0..-4) 

N(4)  =  N(3)  +  6{3)D{3> 

1  v  (0---4)  JV(4)  ~  ^  (4)  ■*t0) 

=  AT(3)  +  i^(3)Z)(0) 

iV(4)  ^/y(4)-Ly(0)  * 


(A.33) 
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The  requirement  for  producing  N^.A)  =  |0)  is  simply  =  ifi$  |-) ;  i.e., 

q0)  =\(-2s4  +\{s,s}  +  \82e" /dx2^ 
r<3>  =-i{k\.£A£,£]'\  +  \{£1  ,£}  +  \{£,de2  Idx}  +  1  ds3 /dx-±d3s/dx3) 


(A.34) 


which  reduces  to  Q'"  =  \a2N\l\ .  This  satisfies  the  pseudounitarity  condition,  and,  it  would  provide  matrix 
representations  for  any  remaining  terms  in  Eq.  (A.33);  however,  all  of  them  are  already  represented.  At 
this  stage,  the  generator  is  diagonalized  through  4th  order  ( 4)  =  |o) );  however,  it  is  still  unbalanced 


Z)t4> 

^(0-4) 


/ 

_ 1_ 

16 

-'■( 

V 

■\Jl  +  2s  j 

-{(£  +  A) 


x\-) 


A.35) 


Balance  is  restored  with  a  block-diagonal  transformation,  Pr"  :  Gi4)  —>  G(4) .  The  form  that  does  this  is 
s  ~  s3 .  (As  with  the  previous  case,  other  options  s  ~  sn  with  n  *  3  lead  to  no  solution.)  With  s  ~  s\  Eqs. 
(62)  and  (63)  yield 


(4)  =  -(4)  ,r  1,7(4) 

w(0---4)  w(3---4)  ^  L*5  Jw(l) 

=  dw  +s 

“(0-4)  "(0-4)tj 


(A.36) 


The  requirement  for  a'4).4)  =  0  is  simply  «‘44’  +[.v|J'4’  =  0 ;  i.e.,  j^[s2,s]  =  [£,w] .  From  the  general  3rd-order 
form  is  =  /F  +  rjss  +  y/ss  +  6'i: ,  we  find  that  6  must  vanish  and  that  t//  =  tj  =  3/1 6  ,  although  %  remains 
undetermined  (because  f^,£-3]  =  0).  Thus,  is  =  and  is  =  xds3/dx  +  ^s2  {s,s} .  However, 

left-right  symmetry  requires  j  =  0  to  eliminate  the  ce  jcx  term.  Thus,  through  4th  order,  the  transformed 
generator  remains  diagonal  ( N^„4)  =  |o) )  and  is  also  balanced: 

d$-4)  =  (-' (■ Vl  +  2e )(o  4)  (A.37) 


The  contributions  that  do  not  belong  to  the  Taylor  series  for  \J\  +  2e  include  the  3rd-order  one  that  persists 
from  the  previous  step  and  a  pair  of  new  4th-order  terms. 


